/*******************************************************************
 * Copyright (C) yel_hb-USTC, 2012                                
 *                                                                  
 * FILE NAME             :  kmeans.cpp                                 
 * PRINCIPAL AUTHOR      :  yel_hb                               
 * MODULE NAME           :  KMEANS                             
 * LANGUAGE              :  C++                                      
 * TARGET ENVIRONMENT    :  windows                                   
 * DATE OF FIRST RELEASE :  2012/12/21                              
 * DESCRIPTION           :  implementaion the methods defined in
 *							kmeans.h
 ******************************************************************/

/*
 * Revision log:
 *
 * Created by yel_hb,2012/12/21
 * Modified by yel_hb,2012/12/23
 * Modified by yel_hb,2012/12/24
 *
 */

#include <vector>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include "dmvector.h"
#include "kmeans.h"
#include "kahansum.h"
#include "utils.h"
using namespace std;

KMeans::KMeans()
{
	k = 1;
	cluster_centers.clear();
	point_ids.clear();
}

KMeans::KMeans ( int k )
{
	this->k = k;
	cluster_centers.resize ( k );
	point_ids.clear();
}

int KMeans::clusterize ( vector<DMVector> vlist )
{
	int itercnt = 0;
	int chcnt = 0x7FFFFFFF;
	int cid, i, j, size;
	vector<double> vb;
	vector<DMVector> vcls;

	size = vlist.size();            //sum of points

	//the criteria of convergence is none point changes
	while ( chcnt != 0 )
	{
		//find the closest clustercenter
		chcnt = 0;
		for ( i = 0; i < size; ++i )
		{
			cid = get_closest_center_id ( vlist[i] );
			if ( get_point_id ( i ) != cid )
			{
				point_ids[i] = cid;
				chcnt++;
			}
		}
    
		//Calculate the initial center
		for ( i = 0; i < k; ++i )
		{
			vcls.clear();
			vb.clear();
			for ( j = 0; j < size; ++j )
			{
				if ( point_ids[j] == i )
				{
						vcls.push_back ( vlist[j] );
				}
			}

			vb = vector_mean ( vcls );
			for ( j = 0; j < vb.size(); ++j )
			{
				cluster_centers[i].set ( j, vb[j] );
			}
		}

		itercnt++;
	}

	return itercnt;
}

int KMeans::get_point_id ( int i )
{
	return point_ids[i];
}

DMVector KMeans::get_cluster_center ( int i )
{
	return cluster_centers[i]; 
}

KMeans::~KMeans()
{
	k = 0;
	cluster_centers.clear();
	point_ids.clear();
}

int KMeans::get_closest_center_id ( const DMVector& x )
{
	return get_closest_center_id ( x, cluster_centers.size() );
}

int KMeans::get_closest_center_id ( const DMVector& x, int size )
{
	int i, mini;
	double mins = 0x7FFFFFFF;
	double len;

	mini = 0;
	for ( i = 0; i < size; ++i )
	{
		DMVector tmp = x;
		tmp.sub ( cluster_centers[i] );
		len = tmp.norm ();
		if ( len < mins )
		{
			mins = len;
			mini = i;
		}
	}
	return mini;
}

void KMeans::init_point_indexes()
{
	vector<int>::iterator ix;
	srand( (unsigned)time( NULL ) );
	for ( ix = point_ids.begin(); ix != point_ids.end(); ++ix )
	{
		(*ix) = rand() % k;
	}
}

void KMeans::init_cluster_center ( vector<DMVector> vlist )
{
	int i, j, size;
	vector<double> vb;
	vector<DMVector> vcls;

	point_ids.resize ( vlist.size() );
	init_point_indexes();
	size = vlist.size();

	for ( i = 0; i < k; ++i )
	{
		vcls.clear();
		vb.clear();
		for ( j = 0; j < size; ++j )
		{
			if ( point_ids[j] == i )
			{
				vcls.push_back ( vlist[j] );
			}
		}

		vb = vector_mean ( vcls );
		for ( j = 0; j < vb.size(); ++j )
		{
			cluster_centers[i].set ( j, vb[j] );
		}
	}
}

void KMeans::kmpp_init ( vector<DMVector> vlist )
{
	int i, j, size, clst;
	double sum;
	DMVector vc;
	vector<double> vd;
	KahanSum kahan;
	srand( (unsigned)time( NULL ) );
	point_ids.resize ( vlist.size() );
	init_point_indexes();

	size = vlist.size();
	vc = vlist[rand() % size];

	for ( i = 0; i < vc.size(); ++i )
	{
		cluster_centers[0].set ( i , vc.get ( i ));
	}

	for ( i = 1; i < k; ++i )
	{
		kahan.reset();
		for ( j = 0; j < size; ++j )
		{
			clst = get_closest_center_id ( vlist[j], i );
			vc = vlist[j];
			vc.sub ( cluster_centers[clst] );
			vd.push_back ( vc.norm() );
			kahan.add ( vc.norm() );
		}
		sum = kahan.get_sum() * rand() / ( RAND_MAX - 1 );

		for ( j = 0; j < vd.size(); ++j )
		{
			if ( sum -= vd[j] > 0 )
				continue;
			else
			{
				vc = vlist[j];
				for ( k = 0; k < vc.size(); ++k )
				{
					cluster_centers[i].set ( k , vc.get ( k ));
				}
			}
		}
	}
}

double KMeans::DI( vector<DMVector> vlist )
{
	vector<double> vmaxi, vmaxj;
	KahanSum khssum;
	DMVector vc;
	double dmax;
	int i, j, size;

	size = vlist.size();
	vmaxj.clear();
	for ( i = 0; i < k; ++i )
	{
		vmaxi.clear();
		for ( j = 0; j < size; ++j )
		{
			if ( get_point_id ( j ) == i )
			{
				vc = vlist[j];
				vc.sub ( cluster_centers[i] );
				vmaxi.push_back ( vc.norm() );
			}
		}
		vmaxj.push_back ( *max_element ( vmaxi.begin(), vmaxi.end() ));
	}
	dmax = *max_element ( vmaxj.begin(), vmaxj.end() );
	
	vmaxj.clear();
	for ( i = 0; i < k; ++i )
	{
		vmaxi.clear();
		for ( j = 0; j < k; ++j )
		{
			if ( i != j )
			{
				vc = cluster_centers[j];
				vc.sub ( cluster_centers[i] );
				vmaxi.push_back ( vc.norm() );
			}
		}
		vmaxj.push_back ( *min_element ( vmaxi.begin(), vmaxi.end() ));
	}

	return ( *min_element( vmaxj.begin(), vmaxj.end() )) / dmax;
}

double KMeans::DBI( vector<DMVector> vlist )
{
	int i, j, s, size, cnti, cntj;
	KahanSum khsmax;
	KahanSum khssumi, khssumj;
	DMVector vc;
	
	khsmax.reset();
	size = vlist.size();
	for ( i = 0; i < k; ++i )
	{
		khssumi.reset();
		cnti = 0;
		for ( s = 0; s < size; ++s )
		{
			if ( get_point_id ( s ) == i )
			{
				vc = vlist[s];
				vc.sub ( cluster_centers[i] );
				khssumi.add ( vc.norm() );
				cnti++;
			}
		}

		for ( j = 0; j < k; ++j )
		{
			if ( i != j )
			{	
				khssumj.reset();
				cntj = 0;
				for ( s = 0; s < size; ++s )
				{
					if ( get_point_id ( s ) == j )
					{
						vc = vlist[s];
						vc.sub ( cluster_centers[j] );
						khssumj.add ( vc.norm() );
						cntj++;
					}
				}
				
				vc = cluster_centers[i];
				vc.sub ( cluster_centers[j] );
				khsmax.add (( khssumi.get_sum() / cnti + khssumj.get_sum() / cntj ) / vc.norm() );
			}
		}
	}

	return khsmax.get_sum() / k;
}

